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1.  Introduction 


1.1  Background 

The  U.S.  Army  Research  Laboratory  (ARL)  in  Aberdeen,  MD,  has  expressed  a  need  for  a 
computer  software  package  that  can  analyze  ion-implantation  data,  extract  diffusion  coefficients 
from  this  data,  and  predict  concentration  profiles.  Their  current  implantation  data  is  used  to 
study  the  effects  of  carbon  diffusion  in  gun-tube  barrels.  A  comprehensive  software  package 
such  as  this  would  save  a  tremendous  amount  of  experimentation  time  and  money  for  ARL  and 
would  therefore  be  very  beneficial. 

1.2  Problem  Statement 

The  objective  of  this  thesis  was  to  develop  a  Fortran  software  package  in  order  to  extract 
diffusion  data  from  concentration  profiles  and  to  predict  future  concentration  profiles.  This 
package  needed  to  be  stand-alone,  user  friendly,  and  have  the  ability  to  interface  with  ARL’s 
existing  computer  code. 

1.3  Approach 

The  approach  to  this  problem  began  with  an  in-depth  literature  survey.  This  allowed  the 
gathering  of  necessary  equations  and  solution  methods  to  solve  the  mathematical  portion  of  the 
thesis.  Using  this  newly  acquired  knowledge,  a  Fortran  code  was  written  in  order  to  satisfy  the 
requirements  of  ARL  and  to  perform  the  necessary  functions.  This  code  was  tested  and 
compared  against  existing  code  and  mathematical  solutions  to  validate  its  effectiveness. 

The  thesis  then  covered  various  experiments  to  demonstrate  the  sensitivity  of  the  code.  The 
purpose  of  this  step  was  to  aid  the  reader  in  visualizing  the  adverse  effects  of  altering  parameters 
such  as  number  of  data  points,  time  increments,  and  induced  noise. 

The  thesis  concluded  with  a  discussion  on  the  limitations  of  the  software  itself  and  possible 
recommendations  that  can  be  made  in  the  future  to  improve  the  software  package. 


2.  Literature  Survey 

2.1  Diffusion  Equations  for  Predicting  Concentration  Profiles 

2.1.1  Binary  Systems 

2. 1 . 1 . 1  Constant  Diffusivity.  In  reality,  most  cases  of  diffusion  are  transient  or  non  steady-state 
ones.  This  applies,  for  example,  to  those  cases  in  which  the  interstitial  concentration  C  varies 
with  time  which  results  in  a  net  accumulation  or  depletion  of  the  diffusing  species.  In  order  to 
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model  this  situation,  it  is  necessary  to  use  the  partial  differential  equation  known  as  the  diffusion 
equation  which  is  given  by  (1 ) 


ac __a(  acN 

dt  dx  j 


(1) 


For  cases  in  which  D  is  independent  of  composition,  or  where  the  range  of  composition  is  small, 
equation  1  reduces  to  (1 ) 


dC  _Dd2C 

dt  dx2 


(2) 


which  is  known  as  Fick’s  second  law.  Figure  1  shows  an  example  of  transient  diffusion  in  a 
binary  system  in  which  the  diffusivity  is  constant  with  respect  to  composition. 


Figure  1.  Transient  diffusion  with  constant  diffusivity. 

If  one  were  to  apply  Fick’s  second  law  to  a  semi-infinite  solid  and  held  the  surface  concentration 
constant,  a  common  error  function  solution  would  be  obtained  (2).  In  order  to  obtain  such  a 
solution,  the  following  assumptions  can  be  made: 

1 .  Before  diffusion,  the  diffusing  solute  atoms  in  the  solid  are  uniformly  distributed  with 
concentration  of  C0. 

2.  The  concentration  at  the  surface,  x  =  0,  is  a  constant  value,  Cs 
These  boundary  conditions  can  be  stated  as  follows: 

For  t  =  0,  C  =  Co  at  0  <  x  <  oo, 


For  t  >  0,  C  =  Cs  at  x  =  0,  and  C  =  Co  at  x  =  oo. 
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If  these  boundary  conditions  are  applied  to  equation  2,  the  following  solution  is  obtained: 


C  -C 
C  -C 

s  o 


(3) 


where  Cx  represents  the  concentration  at  depth  x  after  time  t.  The  expression  erf(x/2/ Dt)  is  the 
Gaussian  error  function,  values  of  which  are  given  in  mathematical  tables  for  various  x/2/ Dt 
values.  Equation  3  demonstrates  the  relationship  between  concentration,  position,  and  time, 
namely,  that  Cx,  being  a  function  of  the  dimensionless  parameter  x// Dt,  may  be  determined  at 
any  time  and  position  if  the  parameters  C0,  Cs,  and  D  are  known. 


Another  solution  to  Fick’s  second  law  is  known  as  the  thin-film  solution.  Taking  an  infinite 
plane  as  the  geometry,  the  total  amount  of  substance  M  diffusing  in  the  cylinder  and  unit  cross 
section  is  given  by  (1 ) 


/•  oo 

M  =  [  Cdx . 

J  —00 


(4) 


After  differentiating  equation  4  and  applying  the  appropriate  derivation,  equation  5  is  given  as 

V) 


exp 


(5) 


Therefore,  this  is  the  solution  which  describes  the  spreading  by  diffusion  of  an  amount  of 
substance  M  deposited  at  time  t  =  0  in  the  plane  x  =  0.  Figure  2  shows  typical  distributions  at  six 
successive  times. 


Figure  2.  Concentration-distance  curves  for  an  instantaneous  plane  source. 
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For  the  thin-film  solution,  we  can  consider  the  solution  for  negative  x  to  be  reflected  in  the  plane 
x  =  0  and  superimposed  on  the  original  distribution  in  the  region  x  >  0.  Since  the  original 
solution  was  symmetrical  about  x  =  0  the  concentration  distribution  for  the  semi-infinite  plane  is 
given  by 


C  = 


M 

V 7lDt 


exp 


^-x2^ 


4  Dt 


(6) 


A  typical  concentration  distribution  for  the  thin  film  solution  is  demonstrated  in  figure  3. 


Figure  3.  Concentration-distance  curves  for  the  thin-film  solution. 


A  third  solution  to  Fick’s  second  law  is  based  on  a  trigonometric  solution.  In  this  solution,  it  is 
assumed  that  the  time  and  spatial  variables  are  separable,  that  is  (7) 

CB=4(x)-T(t),  (7) 


where  ^  and  r  are  functions,  as  yet  unknown,  of  x  and  t  respectively.  Substituting  for  C«  in 
Fick’s  second  law  gives  ( 1 ) 


£ 


dr 

dt 


=  Dbt 


d2% 

dx2 


(8) 


1  dr  1  d2d, 
Dbt  dt  %  dx2 


(9) 
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where  we  now  have  only  total  differentials  of  our  unknown  functions.  Since  one  side  of  the 
equation  depends  only  on  t  and  the  other  only  on  x,  the  equation  can  only  hold  for  all  x  and  t  if 
both  are  equal  to  some  constant,  say  -k  (1 ). 

rj  T 

-  =  -k2DBr.  (10) 

at 

^  =  -k2%.  (11) 

dx 

These  differential  equations  can  now  be  solved  to  give  the  functions  E,  and  z  ( 1 ). 

r  =  exp(-k2DBt)  .  (12) 

£  =  Ax  sin( kx)  +  A2  cos (kx) .  (13) 


The  time-dependent  concentration  function  Cb  is  therefore  a  sinusoidal  composition  fluctuation 
which  decays  exponentially  with  time.  This  is  seen  more  clearly  by  taking  Cb  =  Co  at  x  =  0  at  all 
times,  and  substituting  for  the  wavelength  of  the  sinusoidal  variation  A  : 


CB  (x,  t)  =  C0+  AC  sin 


2  7a 

T~ 


exp 


4  n2DB 
A2 


\ 

t 

y 


=  C0  +  AC  sin 


2  2 7ZK^ 

r  o 

y  A  ) 

exp 

— 

1  j 

(14) 


where 


tT 


A2 

4  n2DB 


(15) 


is  the  relaxation  time,  which  is  the  time  taken  for  a  sinusoidal  variation  of  wavelength  A  to  drop 
to  37%  of  its  original  amplitude. 


2. 1 . 1 .2  Variable  Diffusivity.  When  the  diffusion  coefficient  D  is  a  function  of  concentration  C, 
the  equation  for  one-dimensional  (1-D)  diffusion  is  (3) 


dC__H(DdC' 

dt  dx ; 


(16) 


Differentiation  of  equation  16  yields  (3) 


8C  n  d2C  dD  dC 

—  =  D — -  + - . 

dt  dx  dx  dx 

Second-order-correct  approximations  to  the  partial  derivations  in  equation  17  are  (3): 


(17) 


5 


(18) 


dC 


n+1/2 


dt, 


^ ~in+\  (jn 

At 


2  /^i n+l/ 2 


dzc 


dxt 


s~in+ 1  rss~in+i  .  s~in+\  s~in  _i_ 

L/z+ 1  _  ZC/  |  U/+1~ZC/  '  '“'/-l 


ac 


1+1/2 


M2 

3D-  ,  A’.  - Q-, 
Sx,.  2Ax 

+i 


(A,)2 


(19) 

(20) 


dx, 


s~in+l  s~in+\.  s^in  s^in 
'“'i+l  ~^i- 1  _j_  ^i+l  _  ^i-l 


2Ax 


2Ax 


(21) 


The  superscripts  refer  to  the  time  dimension  and  subscripts  denote  the  space  dimension.  Note 
dD 

that  —  is  evaluated  at  n  rather  than  n+1/2,  because  the  composition-dependent  coefficient  D 

dx 

cannot  be  calculated  at  the  next  time  step  n+l  before  the  concentrations  at  n+l  have  been 
evaluated.  Substituting  the  approximations,  equations  18-21  into  equation  17  yields  (5) 


C':'[-4D-  +D-tl-D’_,]+cr' 


8(Ax)“ 

At 


+  8  D; 


+  C’t'  [-  4£>,"  -  ,  -  D"_t 


ciMd;  -D-t,+Di]+c- 


8(Av)“ 

At 


-8  D: 


+  C+Z>,' +£«-£>■,],  (22) 

which  has  been  arranged  so  that  all  the  concentrations  at  the  current  time  step  (n)  are  on  the  right 
and  the  concentrations  to  be  computed  at  the  next  time  step  (n+l)  are  on  the  left. 

2.1.2  Multicomponent  Systems  With  One  Fast  Diffuser 

The  basis  for  modeling  multicomponent  systems  with  one  fast  diffuser  comes  from  the  following 
equation: 

rs  total  rs 

^ —  =  -(-J  )t0,al  (23) 

dtdx 

in  which  x  represents  distance  and  Jto,al  is  the  total  flux  ( 4 ).  Assuming  the  precipitate  volume 
fraction  is  negligible,  the  multicomponent  flux  can  be  written  as 
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(24) 


n— 1 


•a=-2A 

j= i 


dcj 

dx 


in  which  component  1  is  a  fast  diffuser  (i.e.,  an  interstitial  atom)  and  components  2  through  n-1 
refer  to  slower  moving  alloying  elements  (i.e.,  substitutional  atoms)  (5). 


Equation  24  can  be  rewritten  in  terms  of  an  effective  diffiisivity  described  by  the  equation: 

c)cm 

\eff  OC\ 


J  =  -De 


dx 


(25) 


Joining  equations  24  and  25  yields  (6) 


i  dc" 


eff  _  7=1 


Df  = 


dx 


dc"' 

dx 


(26) 


By  assuming  local  equilibrium  and  no  long-range  diffusion  by  substitutional  atoms,  equation  26 
can  be  simplified  to 


Df  =  D, 


^2,  dc 

7=2  dc 


(27) 


Equation  27  represents  only  one  independent  concentration  variable  with  only  one  degree  of 
freedom. 


Combining  equations  23  and  25  and  then  differentiating  c[ot  gives 

dcT_  =  ^Def/dci_  dcT_ 
dt  dx  1  dc[ot  dx 


(28) 


Equation  28  is  the  one  used  to  predict  concentration  profiles  for  multicomponent  systems  with 
one  fast  diffuser. 


2.2  Finite-Difference  Methods  for  Predicting  Concentration  Profiles 

The  objective  of  a  finite-difference  method  for  solving  an  ordinary  differential  equation  (ODE)  is 
to  transform  a  calculus  problem  into  an  algebra  problem  by: 

1 .  Discretizing  the  continuous  physical  domain  into  a  discrete  finite  difference  grid, 

2.  Approximating  the  exact  derivatives  in  the  initial-value  ODE  by  algebraic  finite  difference 
approximations  (FDAs), 
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3.  Substituting  the  FDAs  into  the  ODE  to  obtain  an  algebraic  finite  difference  equation 
(FDE),  and 

4.  Solving  the  resulting  algebraic  FDE. 

2.2.1  Explicit  Forward-Difference  Method  (Euler  Method) 

Consider  the  general  nonlinear  first-order  ODE  (7): 

y  =/(f»  y)  y(1o)  =  y0-  (29) 

Choose  point  n  as  the  base  point  and  develop  a  finite  difference  approximation  of  equation  29  at 
that  point.  The  finite-difference  grid  is  illustrated  in  figure  4,  where  the  x  symbol  denotes  the 
base  point  for  the  finite  difference  approximation  of  equation  29. 


Figure  4.  Finite  difference  grid  for  the  explicit  Euler  method. 


The  first-order  forward-difference,  finite-difference  approximation  of  y'  is  given  by  (7): 

7|  =-V"+i~-V”  -i-7(r n)At. 

A  t  2 

Substituting  equation  30  into  equation  29  and  evaluating  f[t,y )  at  point  n  yields  (7) 


=  f{tn,y,)= f  n. 


Solving  equation  31  for  yn+1  gives  (7) 

y„+i  =yn+Atfn  +  |v"(r„)Ar  =yn  +Atfn+0(At2).  (32) 

Truncating  the  remainder  tenn,  which  is  o(Ar ),  and  solving  for  yn+1  yields  the  explicit  Euler 
finite  difference  equation  (FDE)  (7): 

=  v„+A(/,o(Ar),  (33) 

where  the  o(a r )  term  is  included  as  a  reminder  of  the  order  of  the  local  truncation  error. 

Several  features  of  equation  33: 

1 .  The  FDE  is  explicit,  since  fn  does  not  depend  on  yn+1 . 
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2.  The  FDE  requires  only  one  known  point.  Hence,  it  is  a  single-point  method. 

3.  The  FDE  requires  only  one  derivative  function  evaluation  (i.e.,  f(t,y))  per  step. 

4.  The  error  in  calculating  yn+]  for  a  single  step,  the  local  truncation  error,  is  o(Ar ) . 

5.  The  global  (i.e.,  total)  error  accumulated  after  N  steps  is  C)(Ar ) . 

The  explicit  Euler  method  only  has  first-order  accuracy  and  it  is  very  unstable,  therefore,  it  is 
impractical  to  use. 

2.2.2  Crank-Nicolson  Method 

The  Crank  Nicolson  (7)  method  is  a  more  widely  used  finite  difference  method  for  solving 
partial  differential  equations  and  is  set  up  using  the  following  grid  (figure  5). 


Figure  5.  The  Crank-Nicolson  method  stencil. 


Crank  and  Nicolson  in  1947  proposed  approximating  the  partial  derivative  ft  at  grid  point 
(i,«+l/2)  by  the  second-order  centered-time  approximation  obtained  by  combining  Taylor  series 
for  /"  and  /”.  Thus,  (7) 


—~n+\  — 22 +1  /  2  — 

ft  =fi  +ft 


1+1/2 


v  *  ) 


+  2f" 


2+1/2 


~~n  — 22+ 1/2  — 

ft  =fi  -ft 


n+U  2 


v  ^  j 


+  2f" 


22+1/ 2 


At 

vTy 

At 


1  —  |  22+1/ 2 

+  -/m 

6 


V  ^  J 


1  _  |„+l/2 

—Jut 

6 


At 

v  2  j 

At 


+  ... 


+  ... 


V  ^  J 


Subtracting  these  two  equations  and  solving  for  f  t 


2+1/2 


gives 


— —22+1  ~  22 

”+1/2_//  -ft  1+  (T\Kf2 
t  At  24  ’ 


(34) 

(35) 


(36) 
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where  t"  <  r  <  t"+l  ^7l  Truncating  the  remainder  term  in  equation  36  yields  the  second-order 
centered-time  approximation  of/,  (7): 


n+ 1/2 


At 


(37) 


The  partial  derivative  f  xx  at  grid  point  (i,n+ 1/2)  is  approximated  by  (7) 


n+ 1/2 


(38) 


The  order  of  the  FDE  obtained  using  equations  37  and  38  is  expected  to  be  0(Af)2  +  0(Ax)“  ,  but 
that  must  be  proven  from  the  MDE.  The  partial  derivative  /  tv  at  time  levels  n  and  n+l  are 
approximated  by  the  second-order  centered-difference  approximation 


V7  +/m 

Ax2 


(39) 


applied  at  time  levels  n  and  n+l,  respectively  (7).  The  resulting  finite-difference  approximation 
of  the  1-D  diffusion  equation  is  (7) 


^*«+ 1  j'n 

At 


1 

-a  — 
2 


f 


n+i  r\  r  n+l  .  r  n+l  rn  r\  rn  .  rn 

i+ 1  ~^Ji  +  J i-\  ,  Ji+l  ~  i  +  //-l 


V 


Ax" 


+  - 


Ax" 


(40) 


Rearranging  equation  40  yields  the  Crank-Nicolson  finite-difference  equation: 

-  d-:; + 2(1 + d)/r‘  -  + 2(1  -  </)/,■ + , 


(41) 


where  d  =  a 


At 

Ax2 


is  the  diffusion  number  (7). 


The  Crank-Nicolson  method  is  unconditionally  stable  and  accurate  on  a  second  order  level.  The 
solution  at  a  given  time  level  can  be  reached  with  much  less  computational  effort  by  taking 
larger  time  steps.  The  time  step  is  limited  only  by  accuracy  requirements. 

2.2.3  Thomas  Algorithm  to  Solve  a  Tridiagonal  System  of  Equations 

When  a  large  system  of  linear  algebraic  equations  has  a  special  pattern,  such  as  a  tridiagonal 
pattern  as  in  the  Crank-Nicolson  equation,  it  is  usually  worthwhile  to  develop  special  methods 
for  that  unique  pattern.  These  methods  are  generally  very  efficient  in  computer  time  and  storage. 
One  algorithm  that  deserves  special  attention  is  the  algorithm  for  tridiagonal  matrices,  often 
referred  to  as  the  Thomas  algorithm. 
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To  derive  the  Thomas  algorithm,  the  Gauss  elimination  procedure  is  applied  to  a  tridiagonal 
matrix  T,  modifying  the  procedure  to  eliminate  all  unnecessary  computations  involving  zeros. 
Consider  the  matrix  equation: 


Tx  =  b,  (42) 

where  T  is  a  tridiagonal  matrix  (7).  Thus,  (7) 


au 

(N 

0 

0 

0 

0 

0 

0 

<N 

a  22 

a  23 

0 

0 

0 

0 

0 

0 

ai2 

«33 

fl34 

0 

0 

0 

0 

0 

0 

°43 

a44 

aA5 

0 

0 

0 

0 

0 

0 

0 

0 

•••  an-\,n-2 

an-l,n~l 

a  , 

n-\,n 

0 

0 

0 

0 

0 

0 

a  i 

n,n- 1 

a 

n,n 

Since  all  the  elements  of  column  1  below  row  2  are  already  zero,  the  only  element  to  be 


eliminated  in  row  2  is  a2l .  Thus,  replace  row  2  by  R2 - 


(  _  A 

<2i 


V  an  J 


Rl .  Row  2  becomes  (7) 


0  a22  - 

fN 

Q  | 

au  a23  0  0. 

..000 

(44) 


Similarly,  only  a32  in  column  2  must  be  eliminated  from  row  3,  only  a43  in  column  3  must  be 

eliminated  from  row  4,  etc.  The  eliminated  element  itself  does  not  need  to  be  calculated.  In  fact, 
storing  the  elimination  multipliers,  em  =  ( a2l  /  an  ),  etc.,  in  place  of  the  eliminated  elements 

allows  this  procedure  to  be  used  as  an  LU  factorization  method.  Only  the  diagonal  element  in 
each  row  is  affected  by  the  elimination.  Elimination  in  rows  2  to  n  is  accomplished  as  follows 
(7): 


a‘-U  (i 


2,...,n) 


(45) 


Thus,  the  elimination  step  involves  only  2n  multiplicative  operations  to  place  T  in  upper 
triangular  form. 


The  elements  of  the  b  vector  are  also  affected  by  the  elimination  process.  The  first  element  b\  is 
unchanged.  The  second  element  becomes  (7) 


b2=b2- 


ra  A 
U21 

V  an  J 


V 


(46) 
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Subsequent  elements  of  the  b  vector  are  changes  in  a  similar  manner.  Processing  the  b  vector 
.....  .  .  .  ...  .  ....  (a,,  . 


requires  only  one  multiplicative  operation,  since  the  elimination  multiplier,  em  =  —  ,  is 

already  calculated.  Thus,  the  total  process  of  elimination,  including  the  operation  on  the  b 
vector,  requires  only  3n  multiplicative  operations. 

The  n  x  n  tridiagonal  matrix  T  can  be  stored  as  an  n  x  3  matrix  A  ’  since  there  is  no  need  to  store 
the  zeros.  The  first  column  of  matrix  A  elements  a, , ,  corresponds  to  the  sub-diagonal  of  matrix 

T,  elements  ai  M  .  The  second  column  of  matrix  A  elements  a\  2 ,  corresponds  to  the  diagonal 

elements  of  matrix  T,  elements  ai . .  The  third  column  of  matrix  A  elements  at  3 ,  corresponds  to 


ai,M  ■ 

The  elements 

a\. 

i  and  an  3  do  not  exist.  Thus,  (7) 

- 

a\,2 

a\,l 

a2,\ 

a2,2 

a2,2 

a2,\ 

a2,2 

a2,3 

(47) 

a'n-u 

an- 1,2 

an- 1,3 

a,ul 

an,2 

- 

When  the  elements  of  column  1  of  matrix  A  ’  are  eliminated,  that  is,  the  elements  at , ,  the 
elements  of  column  2  of  matrix  A  ’  become  (7) 


(/  =  2,  3,  ...,ri) 


The  b  vector  is  modified  as  follows:  (7) 


b,.  =  b,  - 


bi  =bi 


(/  =  2,  3,  ...,ri) 


After  a 2  (i  =  2,  3,  ...,n)  and  b  are  evaluated,  the  back  substitution  step  is  as  follows:  (7) 


(A  ~a[3xi+l) 
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Pivoting  destroys  the  tridiagonality  of  the  system  of  linear  algebraic  equations,  and  thus  cannot 
be  used  with  the  Thomas  algorithm.  Most  large  tridiagonal  systems  that  represent  real  physical 
problems  are  diagonally  dominant,  so  pivoting  is  not  necessary.  The  Thomas  algorithm,  in  a 
format  suitable  for  programming  for  a  computer,  is  summarized  as  follows: 

1 .  Store  the  n  x  n  tridiagonal  matrix  T  in  the  n  x  3  matrix  A  ’.  The  right-side  vector  b  is  an  n  x 
1  column  vector. 

2.  Compute  the  a]  2  terms  from  equations  48  and  49.  Store  the  elimination  multipliers, 
em=  at  ]  /  a)  u ,  in  place  of  ai  t . 

3.  Compute  the  b,  terms  from  equations  50  and  5 1 . 

4.  Solve  for  x,  by  back  substitution  using  equations  52  and  53. 

2.3  Inverse  Methods  Used  to  Calculate  Transport  Properties 


2.3.1  Extracting  Thermal  Conductivity  Values  From  Temperature  Profiles 

Inverse  determination  of  the  thermal  conductivity  from  measured  temperature  profiles  has  been 
the  topic  of  research  by  many  investigators  (3).  Most  of  these  studies  assume  that  the  thermal 
conductivity  is  only  a  function  of  the  spatial  coordinate.  However,  thermal  conductivities  are 
temperature-dependent  quantities  in  most  practical  engineering  applications.  Yeung  developed  a 
second-order  finite-difference  procedure  for  the  inverse  detennination  of  the  thermal 
conductivity  in  a  one-dimensional  heat  conduction  domain.  In  this  case,  the  thermal 
conductivity  of  the  material  is  reconstructed  by  using  the  available  temperature  data  at  discrete 
grid  points.  The  numerical  procedure  is  validated  by  comparing  it  to  known  examples.  It  is 
proven  that,  using  this  technique,  a  priori  knowledge  of  the  functional  form  for  thermal 
conductivity  is  not  required. 


2.3.2  Extracting  Diffusion  Coefficients  From  Concentration  Profiles 

Since  the  governing  equations  for  heat  conduction  and  diffusion  are  similar,  it  is  only  natural  to 
use  the  same  procedure  to  investigate  the  diffusion  coefficient  in  a  concentration  profile. 


As  stated  earlier,  in  a  1  -D  formulation  with  the  diffusing  substance  moving  in  the  direction 
nonnal  to  a  sheet  of  thickness  2a,  the  diffusion  equation  can  be  written  as 


8C_ 

dt 


dC^ 

dx  , 


(0<x<a,  t>0), 


(54) 


where  C  is  the  concentration  of  the  diffusing  substance,  t  is  the  time,  D  is  the  diffusion 
coefficient,  and  x  is  the  distance  coordinate  measured  from  the  center  of  the  sheet  (3). 

Let  the  initial  condition  be  (3) 
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C  =  Co(0<x<a,  t  =  0), 


(55) 


where  Co  is  a  constant  concentration  in  the  medium,  and  let  the  boundary  conditions  be 


—  =  0  (x  =  0,  t  >  0)  (56) 

dx 

D—  =  S(Ce-C)(x*a,  t>0),  (57) 

ox 

where  S  is  the  surface  emission  coefficient  and  Ce  is  the  equilibrium  concentration  (5). 

The  first  step  in  the  inverse  method  is  to  present  a  finite-difference  procedure  for  the  calculation 
of  the  diffusion  coefficient  at  discrete  grid  points.  Let  half  of  the  medium  thickness,  a,  be 
discretized  with  mesh  width  Ax  in  distance  (thickness  direction)  and  At  in  the  time  direction  with 
grid  points  xj  =  j  ■  Ax  (where  j  =  0,  1 , . . . ,n)  and  tt  =  i  ■  At  (where  i  =  0,  1,2...).  The  present 

procedure  will  assume  that  C(x,t)  is  known  at  grid  points  (xj,  t,).  Equation  54  can  then  be 
discretized  as  follows: 


1 .  At  the  surface  grid  point  with  j  =  0  and  i  >  0: 

Applying  forward  difference  to  the  time  derivative  of  equation  54,  we  have  (5) 

r  dCy' 


V  St  j0 


s-ii+l  /~il 

At 


(58) 


Applying  the  central  difference  to  the  distance  derivative,  we  obtain 


"a 

<dA) 

dx 

(  n‘  +  n‘  C‘  -  C‘  C‘  -  C‘  } 
yA  T  o  Li  ^ o  £)  a  i  '-o 

v  2  Ax  0  Ax  y 

Ax 
~2 


where  the  following  has  been  set  in  equation  57: 


dx 


D0a 


J  o 


zx  i  s~i  i 
^1  ~C0 

Ax 


(59) 


(60) 


by  introducing  an  appropriate  constant  a  to  compensate  the  use  of  forward  difference  in  the 
equation,  which  involves  different  errors  than  central  difference  (5).  This  also  pennits  the 
avoidance  of  using  the  unknown  surface  emission  coefficient  S. 

Equating  equations  58  and  59  gives  (5) 

%L  (C“  -c;)=o;(2a-hc; -c;)+d;(c;-c;).  («i) 
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2.  At  an  internal  grid  point  with  0  <  j  <  n  and  i  >  0: 
Here  we  have  (3) 

facY  _  cf  -c) 

v  dt  Jj  At 

and 


8_ 

dx 


VC' 


V 


dx 


D‘+l  +D‘ 


q+1  -c; 

Ax 


D'j+D'j-  i 


C'  -C‘ 

°J-1 


Ar 


Ar 


Equating  equations  62  and  63  yields  (3) 


(62) 


(63) 


2(Ax)2 

At 


(c'“-C‘) 


=  o'.,(c;.,-c')+d;(c;u 


-2c;  + 


(64) 


3.  At  the  center  grid  point  with  j  =  n  and  i  >  0: 

Due  to  symmetry,  we  can  set  C'  ,  =  C'+],  D'_,  =  D'+l ,  and  j  =  n  in  equation  64  to  obtain  (3) 


(Ay)2 

At 


(cr-c:)=DY1(cY1-c:)+D;(cY1 


(65) 


The  next  step  in  the  inverse  method  occurs  if  C\x,t  j  and  C\x,t  +  At  J  are  known  at  evenly  spaced 
grid  points  where  t  is  the  specified  time  and  At  is  the  time  increment,  and  we  are  interested  in 
finding  the  diffusion  coefficient  values  at  the  grid  points.  From  equations  63,  64,  and  65,  we  can 
create  the  following  system  of  linear  equations: 


Ad  =  b, 


(66) 


where  A  is  an  (n+1)  x  (n+1)  matrix  and  d  and  b  are  (n+1)  vectors  (3).  A,  d,  and  b  are  subscripted 
from  0  to  n  as  shown  by  the  following: 
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a0,0  a0,l 

CL  i  q  Cl  |  i  Cl  i  2 


A  = 


an-\,n-2  an-\,n-\  a  n-\ ,n 


"a>~ 

K 

d  = 

b  = 

A,\ 

Pn_ 

(67) 


The  elements  of  d  are  the  unknown  diffusion  coefficient  values  at  the  grid  points,  and  the 
elements  of  A  and  b  are  expressed  as  follows: 

1 .  At  the  surface  grid  point  with  x  =  xo  and  t  =  t :  (3) 


a00  =(2«-l)[c(x0,t)-c(x1,t)]. 

(68) 

o01  =  c{xx,t)-c{xQ,t). 

(69) 

b0  =  ^  ^  [c(x0,t  +  At)  c{xG,t) 

(70) 

At 

At  an  internal  grid  point  with  x  =  xj  (0<j<n)  and  t  =  t :  (3) 

ajj-i  =c{xj_x,t)-c{xj,t). 

(71) 

ajj  =  c(xj+l ,  t)-  2c(xj ,  t)+  c{x  , 

t).  (72) 

aJJ+l  =c[xj+l,t)-c{xj,t). 

(73) 
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b  2  ^  \c(x ,,t  +  At)  c(x , , f )1 . 

j  At  1  x  J  ’  v  1  'J 

(74) 

At  the  center  grid  point  with  x  =  xn  and  t  =  t :  [3] 

an,n-l=C(Xn-ld)-C(xn,t). 

(75) 

an,„  =c{xn_x,t)-c{xn,t). 

(76) 

bn  =  '  J  [c(xn,~t  +  At)  c{xn,t\ 

(77) 

This  system  consists  of  a  tridiagonal  system  of  linear  algebraic  equations.  The  solution  vector  d 
is  the  diffusion  coefficient  vector.  This  system  can  be  solved  using  the  Thomas  algorithm  as 
previously  mentioned. 

3.  Fortran  Software  Development 


3.1  Concentration  Profile  Predictor 
Constant  Diffusivity 

This  program  was  written  to  predict  future  concentration  profiles  from  an  initial  concentration 
profile  and  constant  diffusivity  value.  The  program  begins  by  asking  the  user  for  the  initial 
concentration  profile  file  in  .txt  format.  After  storing  this  array,  the  program  requests  the 
constant  diffusivity  value  from  the  user.  The  final  request  from  the  software  is  the  time  at  which 
the  user  would  like  the  concentration  profile  to  be  predicted.  The  software’s  output  is  both  on 
screen  and  in  a  .txt  file  located  in  the  same  location  as  the  initial  concentration  profile.  The  user 
can  then  easily  import  this  file  into  Excel  to  see  the  new  concentration  profile.  The  user  can  also 
edit  the  position  and  time  increments  within  the  code  itself  to  tailor  the  output  to  their  liking. 

As  discussed  in  section  2. 1.1.1,  one  solution  to  Fick’s  second  law  is  the  error  function  solution. 
One  way  to  test  the  validity  of  this  program  is  to  compare  it  against  this  known  solution, 
demonstrated  by 


C  -C 

x  o 

C  -C 

s  o 


(78) 


This  is  accomplished  by  using  an  initial  concentration  profile  based  on  equation  78  and  setting  a 
constant  diffusivity.  This  diffusivity  value  is  used  both  in  equation  78  and  in  the  program  itself. 
Equal  times  were  chosen  as  well.  The  result  of  this  comparison  is  shown  in  figure  6. 
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Figure  6.  Graph  showing  the  error  function  solution  compared  with  the  Fortran  code  solution. 


One  can  see  that  the  results  match  to  four  significant  figures,  thereby  verifying  the  correct 
operation  of  the  Fortran  program. 


Another  solution  known  as  the  thin-film  solution  was  discussed  in  section  2. 1 . 1 . 1  as  well.  This 
research  compares  this  solution  to  the  Fortran  program  as  well  to  further  validate  its  functional 
use.  The  thin-film  solution  was  given  in  equation  6,  which  is  shown  again  here: 


C  = 


M 

sj  7zDt 


exp 


4  Dt 


(79) 


In  order  to  compare  the  program  to  the  thin-film  solution,  an  initial  concentration  profile  based 
on  equation  79  was  used  as  well  as  a  predetermined  constant  diffusivity.  This  diffusivity  value 
was  used  both  in  equation  79  and  in  the  program  itself.  Equal  time  increments  were  chosen  as 
well.  The  result  of  this  comparison  is  shown  in  figure  7. 


One  can  see  that  the  results  match  to  four  significant  figures,  thereby  verifying  the  correct 
operation  of  the  Fortran  program. 


A  third  solution  of  Fick’s  law  is  based  on  a  trigonometric  solution  and  is  mentioned  in  section 
2. 1 . 1 . 1 .  This  solution  is  based  on  the  following: 


Cg  (x,  t)  =  CQ  +  AC  sin 


f  ItoC 

f  t) 

y  A  ) 

exp 

— 

1  (r  J 

(80) 
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Figure  7.  Graph  showing  the  thin-film  solution  compared  with  the  Fortran  program. 

In  order  to  compare  the  program  to  the  trigonometric  solution,  an  initial  concentration  profile 
based  on  equation  80  was  used  as  well  as  a  predetermined  constant  diffusivity.  This  diffusivity 
value  was  used  both  in  equation  80  and  in  the  program  itself.  Equal  time  increments  were 
chosen  as  well.  The  result  of  this  comparison  is  shown  below  in  figure  8. 

One  can  see  that  the  results  match  to  four  significant  figures,  thereby  verifying  the  correct 
operation  of  the  Fortran  program. 

3.2  Diffusivity  Extractor 
Constant  Diffusivity 

This  program  was  written  to  extract  diffusivity  values  from  two  concentration  profiles.  This 
matrix  of  diffusivity  values  can  then  be  used  to  predict  future  concentration  profiles  with  respect 
to  both  time  and  temperature.  The  program  begins  by  asking  the  user  for  the  two  concentration 
profiles  file  in  .txt  format.  The  array  size,  time  step,  and  time  duration  are  then  inputted.  The 
software’s  output  is  both  on  screen  and  in  a  .txt  file  located  in  the  same  location  as  the 
concentration  profiles.  The  user  can  then  easily  import  this  file  into  excel  to  see  the  graph  of  the 
diffusivity  matrix. 
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Figure  8.  Graph  showing  the  trigonometric  approximation  compared  with  the  Fortran  program. 

In  order  to  test  this  portion  of  the  software,  the  thin-film  solution  was  used  again.  A  time  of 
0.005  s  was  inputted  into  equation  79  in  order  to  generate  one  concentration  profile.  A  second 
profile  was  generated  using  a  time  of  0.00503  s  demonstrating  a  time  step  of  0.00003  s.  In  this 
case,  the  diffusivity  was  set  as  a  constant  value  of  1  to  generate  both  of  these  curves.  In  order  for 
the  program  to  be  operational,  it  would  need  to  extract  a  diffusivity  matrix  with  the  value  1  in 
each  location.  The  resulting  extraction  is  shown  in  graphical  form  in  figure  9. 

As  one  can  see,  the  diffusivity  values  oscillate  at  first  and  finally  converge  to  a  value  of  1 .06, 
which  leaves  a  5%  error  since  1  is  the  true  value.  Since  this  is  within  the  acceptable  tolerance, 
this  portion  of  the  program  is  deemed  operational. 


4.  Software  Limitations 


The  Fortran  programs  written  for  this  research  project  depend  solely  on  input  and  output  external 
files.  These  files  are  currently  in  .txt  format  which  rely  on  exact  formatting  and  data  placement 
specifications.  Slight  alterations  in  either  of  these  variables  will  inherently  affect  the  operation 
of  the  main  program.  If  these  files  are  converted  to  spreadsheet  files  (i.e.,  Excel),  data 
manipulation  and  representation  will  become  more  efficient.  The  programmer  will  need  to  write 
an  SQL  subroutine  to  allow  standard  Fortran  output  to  be  inserted  into  a  spreadsheet. 
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Figure  9.  Diffusivity  matrix  plotted  for  constant  diffusivity  of  1. 

Throughout  this  research,  many  different  programs  were  written  to  solve  different  aspects  of  the 
problem  such  as  binaryeffd,  inverse-rewrite  ,  and  ternary  separate  and  can  be  seen  in 
appendices  A,  B,  and  C  respectively. 

4.1  Binaryeffd 

This  program  was  written  to  act  as  a  concentration  profile  predictor  for  use  with  constant 
diffusivity.  It  uses  the  condensed  fonn  of  Fick’s  second  law  when  ‘D’  is  not  a  function  of 
concentration.  The  program  uses  the  Crank-Nicholson  finite-difference  method  to  iterate  to 
future  unknown  time  steps  to  predict  concentration  profiles  for  any  given  initial  binary  diffusion 
data  set. 

While  this  program  runs  successfully  using  known  diffusion  examples,  the  largest  error  seems  to 
come  from  boundary  condition  determination.  The  program  began  with  Dirichlet  boundary 
conditions  in  which  specify  the  value  of  the  function  at  the  surface  and  the  finite  difference  only 
takes  place  between  them  (i.e.,  i  =  2. .  .n-1). 

T=f(r,t).  (81) 

This  scenario  did  not  work  with  Fortran  so  the  boundary  conditions  were  changed  to  Neumann 
boundary  conditions  which  specify  the  normal  derivative  of  the  function  on  the  surface. 

f^=n  VT  =  f(r,t).  (82) 

on 
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The  Neumann  scenario  proved  effective  when  comparing  the  programs’  output  to  known 
solutions.  This  program  successfully  uses  the  Thomas  algorithm  to  solve  for  the  concentration 
values  at  future  unknown  time  steps  using  the  previous  known  concentration  values. 


4.2  Inverse  rewrite 


This  program  is  designed  to  act  as  a  diffusivity  extractor  which  would  extract  the  composition- 
dependent  interdiffusion  coefficients  from  the  concentration  profiles  in  a  single  diffusion  couple 
The  procedure  is  based  on  the  minimization  of  the  difference  between  the  profiles  calculated  by 
a  finite  difference  scheme  and  the  experimental  profiles  given  by  ARL.  This  program  works 
flawlessly  when  modeled  after  the  thin-film  solution  as  seen  in  figure  10. 


Figure  10.  Thin  film  approximation  using  program. 

This  program  does  not  perform  as  well  when  using  actual  data  from  ARL’s  test  matrix. 
Unfortunately,  due  to  the  variation  in  the  input  data,  the  output  does  not  converge  on  a  specific 
value  as  seen  in  figure  1 1 . 
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Figure  11.  Program  output  using  ARL  input  data. 

Some  effort  was  made  to  computationally  induce  noise  to  simulate  experimental  concentration 
scatter.  A  Gaussian  noise  with  a  standard  deviation  of  0.65%  was  applied  to  the  raw  data  of 
ARL.  Unfortunately,  this  caused  the  program’s  output  to  fluctuate  even  more.  The  inter- 
diffusion  coefficient  values  were  meant  to  be  obtained  along  the  entire  diffusion  path,  instead  of 
only  at  the  intersection  point  of  independent  paths  (Boltzmann-Matano;  BZMA)  (5)  or  instead  of 
mean  coefficient  values  (Krishtal;  KMAZ)  (P). 

The  researcher  believes  the  error  in  this  program  focuses  on  the  initial  values  that  are  introduced. 
It  is  assumed  that  a  genetic  algorithm  would  need  to  be  developed  to  apply  to  the  first  iterative 
step.  This  algorithm  would  aid  in  the  location  of  a  true  minimum  and  not  a  local  one  as  seen  in 
the  current  program.  It  is  also  alleged  that  the  stopping  criterion  for  this  program  was  not 
developed  properly.  This  criterion  should  allow  a  point  at  when  it  is  reached,  a  test  would  be 
performed  on  the  different  terms  to  ensure  they  are  of  the  same  order  of  magnitude.  In  order  to 
maintain  this  criterion  constant  during  the  while  profile  treatment,  and  in  order  to  limit 
calculation  time,  a  variable  increment  would  need  to  be  implemented.  This  increment  could  be 
used  as  follows:  if  the  relative  variation  of  the  concentration  exceeds  10%,  the  increment  is 
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decreased  so  that  the  relative  variation  is  no  more  than  1%.  In  that  case,  for  the  next  iteration,  it 
is  necessary  to  employ  the  previously  discussed  genetic  algorithm  to  determine  the  interdiffusion 
coefficients.  This  may  explain  the  wide  fluctuations  seen  in  figure  1 1 . 


The  researcher  hypothesizes  that  a  smoothing  subcode  would  need  to  be  applied  on  the  raw  ARL 
date  before  using  the  Fortran  program  on  it.  An  appropriate  smoothing  function  is  described  in 
equation  83. 


c,  W=4,  +  Z - 

1  +  (l-£,,  ■'»!,)'  exp| 


'  x~Pij ^ 
r>J  J 


+  B„  •  exp 


•( x-Py)' 


(83) 


In  order  to  use  the  smoothing  procedure  efficiently,  an  adequate  number  of  functions  (q)  must  be 
used.  In  this  case,  five  should  be  sufficient.  The  average  relative  error  on  the  concentration 
associated  with  the  smoothing  procedure  amounts  to  0.02%,  and  the  maximum  relative  error, 
generally  observed  at  extrema,  never  exceeds  0.5%. 

Inverse  solutions  are  known  to  be  sensitive  to  changes  in  input  data  resulting  from  measurement 
and  modeling  errors.  Hence,  they  may  not  be  unique.  Mathematically,  the  inverse  problems 
belong  to  the  class  of  ill-posed  or  ill-conditioned  problems;  that  is,  their  solutions  do  not  satisfy 
the  general  requirements  of  existence,  uniqueness,  and  stability  under  small  changes  to  the  input 
data  (10). 

4.3  Ternaryseparate 

This  program  was  written  to  scale  the  previous  programs  into  ternary  and  multicomponent 
diffusion  with  constant  and  variable  diffusivity.  In  the  ternary  program,  Fortran  would  not 
recognize  the  second  element  array  no  matter  how  it  was  represented.  A  third  element  is  not 
needed  in  the  program  as  the  third  element  is  always  determined  by  the  balance  of  the  other  two 
to  total  100%.  It  is  here  that  a  two-dimensional  (2-D)  array  written  in  C++  would  be  easier  to 
code  and  debug.  Fortran’s  limitations  make  it  difficult  for  the  program  to  access  a  2-D  matrix 
and  perfonn  the  necessary  calculations  while  keeping  track  of  each  data  point  within  the  matrix. 
When  forming  the  three-dimensional  array  known  as  ‘C’,  Fortran  continually  disallowed  a  non¬ 
integer  in  the  ‘countreal’  slot.  This  variable  began  as  an  integer  but  somehow  was  converted  to  a 
non-integer  within  the  program  itself. 

Versions  of  the  Crank-Nicholson  method  and  the  Thomas  algorithm  were  again  used  in  this 
program  with  adjustments  made  to  account  for  the  additional  variables  used  in  multicomponent 
diffusion.  The  finite-difference  nomenclature  is  the  same  as  that  used  for  the  binary  diffusion 
case.  This  program  had  difficulty  interpreting  the  bi-tridiagonal  matrix  that  is  generated  from 
attempting  to  solve  multicomponent  diffusion.  Only  a  single  tridiagonal  matrix  is  generated 
when  solving  a  binary  diffusion  problem. 
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5.  Conclusions 


In  conclusion,  the  original  objective  of  this  project  was  not  met  to  its  completion.  Several  factors 
contributed  to  this  shortcoming,  but  the  primary  obstacle  was  the  correlation  between  the 
software  and  the  input  data.  While  the  software  ran  successfully  with  many  different  known 
solutions,  it  did  not  perform  well  using  actual  concentration  profile  data  from  ARL.  Most  likely, 
this  is  due  to  the  limited  amount  of  species  data,  the  accuracy  of  the  data  itself,  and  the  spacing 
between  each  data  point.  As  with  all  software,  the  closer  the  data  points  and  the  less  intense  the 
noise  is,  the  more  accurate  the  solution  will  be.  This  is  especially  true  when  using  a  method  such 
as  the  finite  difference  method,  which  is  based  solely  on  iteration. 

In  order  to  make  this  program  successful,  careful  analysis  of  the  input  data  must  be  made.  All 
efforts  should  be  taken  to  obtain  more  accurate,  smoother  input  data  which  will  allow  the 
software  to  run  with  fewer  obstacles  and,  in  turn,  produce  cleaner  output  data.  Once  the  input 
data  is  appropriate,  the  programmer  should  return  to  “fine  tune”  the  individual  software 
programs  to  allow  them  to  work  with  the  new  data.  This  may  require  using  a  filtering  subroutine 
in  order  to  accept  only  worthy  data  from  the  input  stream.  There  should  be  a  minimum  number 
of  data  points  that  the  program  can  run  on  and  still  produce  acceptable  results.  It  is  imperative 
that  this  number  is  determined  so  the  program  can  be  modified  to  run  on  a  sufficient  amount  of 
data  points.  This  part  of  the  code  would  terminate  the  program  if  there  were  fewer  data  points 
than  the  program  needed. 

Finally,  the  concentration  profile  predictor  should  also  be  expanded  to  account  for  variable 
diffusivity. 
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Appendix  A.  Binaryeffd 


Binaryeffd 


module  setupinfo 

real  : :  time,  dt,  dx,  t 
real  ::  alpha,  xpos,  cone 
integer ::  i 
real : :  j 

data  dt,dx,time  /0. 00005, 0.01, .01/ 

end  module  setup  info 

program  binary  effd 

use  setup  info 
implicit  none 
interface 

subroutine  fill_a(C,a,nx,effd) 
integer ::  nx 

real,  dimension(0:nx)  ::  C 
real,  dimension(0:nx,4)  ::  a 
real,  dimension(0:nx)  ::  effd 
end  subroutine  fill  a 

subroutine  tridiag(C,a,nx) 
integer ::  nx 

real,  dimension(0:nx)  ::  C 
real,  dimension(0:nx,4)  ::  a 
end  subroutine  tridiag 

end  interface 


integer,  parameter  ::  NSEG  =  40 


real,  dimension  (0:NSEG)  ::  effd  !This  is  the  efective  D  array 
real,  dimension  (0:NSEG,4)  ::  a 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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real,  dimension  (0:NSEG)  ::  C  [Initial  condition 
real,  dimension  (0:NSEG)  ::  xposa,conca 
integer ::  count 

[character  (len=l)  ::  tab  =  char(9) 

!***  Enter  data  and  info 

open  (unit=10,  file-conc.txt',  status='old')  [Opening  the  initial  concentration  profile 
open  (unit=13,  file- ds.txt',  status='old')  [Opening  the  effective  D  array  file 
open  (unit=100,  file-binaryoutputtrig.txt',  status='unknown') 

!dx=l/NSEG 
[print*,  'dx  is’,dx 
alpha=dt/(dx*dx) 

!write(*,*)dx,dt,'  alpha  is',  alpha 

count=-l; 

5  read  (10,*,END=15)  xpos,conc 
count=count+l; 
xposa(count)=xpos; 
conca(count)=conc; 
go  to  5 

15  if  (count.EQ.O)  then 

print*,  'No  data  in  file’ 

else 
end  if 

C=conca; 

[print*  ,C 


[Input  effective  D’s  into  ’effd'  array 

count=-l; 

6  read  ( 1 3 ,  *  ,END=  1 6)  effd 

count=count+l; 
go  to  6 

16  if  (count.EQ.O)  then 

[print*,  ’No  data  in  file’ 

else 

end  if 

[print*,  effd(O) 


!***Crank  Nicholson  Method 
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t=  100.0 
count=0 


do  j=0,time,0.005 

!  if  (t>time)  exit 

!count=count+l 

!t=count*dt 

call  fill_a(C,a,NSEG,effd) 

call  tridiag(C,a,NSEG)  lUpdate  C(i)  to  new  time  step 


end  do 

print*,  C 
write  (100,99)  C; 

99  format  (IX,  FI 0.4); 

close  (unit=100) 

end  program  binary  effd 


i ubroutine  F ill  ^  h®*!**!**!**!®*!*'!®*!**!* 


subroutine  fill_a(C,a,nx,effd) 
use  setup  info 
integer ::  nx 

real,  dimension(0:nx)  ::  C 
real,  dimension(0:nx,4)  ::  a 
real,  dimension(0:nx)  ::  effd 


do  i=l,nx-l 

a(i,l)  =  -alpha/2.  *effd(i) 
a(i,2)  =  1.  +  (alpha*  effd(i)) 
a(i,3)  =  a(i,l) 

a(i,4)  =  C(i)*(l.-(alpha*effd(i)))  +  ((C(i-l)+C(i+l))*alpha/2.*effd(i)) 

end  do 

Iprint*,  a(:,4); 

!  This  is  for  i=0 
a(0,l)  =  0.*effd(0) 
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a(0,2)  =  1.  +  (alpha  *effd(0)) 
a(0,3)  =  -alpha  *effd(3) 

a(0,4)  =  C(0)*(l.-(alpha*effd(0)))  +  (C(l)*alpha*effd(0)) 

!  This  is  for  i=NSEG 
a(nx,l)  =  -alpha  *effd(nx) 
a(nx,2)  =  1.  +  alpha*  effd(nx) 
a(nx,3)  =  0.*effd(nx) 

a(nx,4)  =  C(nx)*(l.-(alpha*effd(nx)))  +  (C(nx-l)*alpha*effd(nx)) 
Iprint*,  a(nx,4); 
end  subroutine  fill  a 


|  I' matrix 

subroutine  tridiag(C,a,nx) 

real,  dimension(0:nx)  ::  C 
real,  dimension(0:nx,4)  ::  a 
integer ::  nx 
real ::  denom 
integer : : j 

C(0)=a(0,4)/a(0,2) 

a(0,3)=a(0,3)/a(0,2) 

do  j=l,nx 

denom  =  a(j,2)  -  a(j,l)*a(j-l,3) 

CO)  =  (a0\4)  -  a(j ,  1  )*C(j - 1 ))/ denom 
a(j,3)  =  a(j,3)/denom 

end  do 

!a(nx,3)  is  0 
do  j=nx-l,0,-l 

CG)  =  CG)-aG,3)*CG+l) 

end  do 

end  subroutine  tridiag 
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Appendix  B.  Inverserewrite 


Inverse  rewrite 


module  setupinfo 

real  : :  time,  dt,  dx,  t 

real  ::  alpha,  xpos,  cone,  xpos2,  conc2 

integer ::  i 

real : :  j 

data  dt,dx,time  /0. 00005, 0.01, 10/ 


end  module  setup  info 


program  inverse_rewrite 

use  setup  info 
implicit  none 
interface 

subroutine  fill_a(a,b,c,r,nx,yimp,yl50) 
integer ::  nx 

!real,  dimension(0:nx)  ::  C 
real,  dimension(0:nx)  ::  a,b,c,r,h 
!real,  dimension(0:nx,0:nx)  ::  a 
real,  dimension(0:nx)  ::  yimp,yl50 
end  subroutine  fill  a 

subroutine  tridiag(nx,a,b,c,r,h) 
integer ::  nx 

real,  dimension(0:nx)  ::  h,a,b,c,r 
!real,  dimension(0:nx,0:nx)  ::  a 
end  subroutine  tridiag 

end  interface 


integer,  parameter  ::  NSEG  =  67 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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!real,  dimension  (0:NSEG,0:NSEG)  ::  a 
real,  dimension  (0:NSEG)  ::  h,a,b,c,r 
real,  dimension  (0:NSEG)  ::  xposa,yimp,xpos2a,yl50 
integer ::  count 

!  character  (len=l)  ::  tab  =  char(9) 


!***Enter  data  and  info 

open  (unit=10,  file-deltal.txt',  status='old') 
open  (unit=13,  file='delta2.txt',  status='old') 
open  (unit=100,  file='output.txt',  status='unknown') 

!dx=l/NSEG 
Iprint*,  'nx  is',NSEG 
alpha=(dx*dx)/dt 
!write(*,*)dx,dt,’  alpha  is',  alpha 

count=-l; 

5  read  (10,*,END=15)  cone 
count=count+l; 

!  xposa(count)=xpos; 
yimp(count)=conc; 
go  to  5 

15  if  (count.EQ.O)  then 

print*,  'No  data  in  file’ 

else 
end  if 

!  print*,  yimp 

count=-l; 

6  read  (13,*,END=16)  conc2 
count=count+l; 

!xpos2a(count)=xpos2; 

y  1 50(count)=conc2; 
go  to  6 

16  if  (count.EQ.O)  then 

! print*,  'No  data  in  file' 
else 
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end  if 


Iprint*,  y!50 


!***Crank  Nicholson  Method 

t=  100.0 
count=0 

!do  j=0, time, 0.5 

!  if  (t>time)  exit 

!count=count+l 

!t=count*dt 

call  fill_a(a,b,c,r,NSEG,yimp,yl50) 
call  tridiag(NSEG,a,b,c,r,h) 

Iprint*, 'hnx  is',h(NSEG) 

lend  do 

print*,  h 
write  (100,99)  h; 

99  format  (IX,  FI 0.4); 

close  (unit=100) 

end  program  inverse_rewrite 


|  ^ ubrOlitlllC  F ill  3, 

subroutine  fdl_a(a,b,c,r,nx,yimp,yl50) 
use  setup  info 
integer ::  nx 

real,  dimension(0:nx)  ::  a,b,c,r 
Ireal,  dimension(0:nx,0:nx)  ::  a 
real,  dimension(0:nx)  ::  yimp,yl50 

i 'pjug  jg  for  2 — 0 


!a(l)=0!!l 

b(0)  =  2.0*yimp(0)-3.0*yimp(l)+yimp(2) 
Iprint*, b(0) 
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c(0)  =  -yimp(0)+yimp(l) 
!print*,c(0) 

r(0)  =  alpha*(yl50(0)-yimp(0)) 
!print*,'y  150(0)  is',y  1 50(0) 

! print*, 'alpha  is', alpha 
!print*,r(0) 


!*!i5***This  is  for  0<i<NSEG******** 
do  i=l,nx-l 

a(i)  =  yimp(i-l)-yimp(i+l) 

b(i)  =  4*(yimp(i+l)-2*yimp(i)+yimp(i-l)) 

c(i)  =  yimp(i+l)-yimp(i-l) 

r(i)  =  4*alpha*(yl50(i)-yimp(i)) 

end  do 

!print*,'yimp(0)  is',yimp(0) 

! print*, 'yimp(l)  is',yimp(l) 

!print*,'yimp(2)  is',yimp(2) 

!  print*, a(l); 

!  print*  ,b(l); 


I ^ ^ ^ ^ ^ ^  for 
!a(nx,l)  =  0 

b(nx)  =  2*(yimp(nx-l)-yimp(nx)) 
!a(nx,3)  =  0 

r(nx)  =  alpha*(yl50(nx)-yimp(nx)) 
Iprint*,  a(nx,4); 
end  subroutine  fdl  a 


|  XX13.tXlX 

subroutine  tridiag(nx,a,b,c,r,h) 
!real,  dimension  (0:nx,0:nx)  ::  a 
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real,  dimension(0:nx)  ::  a,b,c,r,h 

real  bet 

real  gam(lOO) 

integer j 

if(b(0).eq.0)pause  'tridiag: rewrite  equations' 


bet=b(0) 

Iprint*,  'bet  is’,bet 

h(0)=r(0)/bet 

h(nx)=1.0 

Iprint*, ’h(O)  is’,h(0) 

!  ^Decomposition  and  forward  substitution 
do  2  j=l,nx 

gam(j)=c(j- 1  )/bet 

!print*,gam(l) 

bet=b(j)-a(j)*gam(j) 

!  print*, bet 

if  (bet.eq.O)pause  'tridiag  failed’  !  Algorithm  fails 
hG)=(rG)-aG)*h(j-l))/bet 

2  continue 

!  *  *Backsubstitution*  * 

do  3  j=nx-l,0,-l 

h(j)=h(j)-gam(j+l  )*h(j+l) 

3  continue 
return 

end  subroutine  tridiag 
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Intentionally  left  blank. 
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Appendix  C.  Ternary  separate 


T  ernaryseparate 


module  setupinfo 

real  : :  time,  dt,  dx,  t 

!Xpos  is  the  array  of  positions  in  the  initial  profile 
!Conc  is  the  array  of  concentrations  in  the  initial  profile 
real  ::  alpha,  xposl,  concl,  xpos2,  conc2,dtt 
integer ::  i,k,z 
real : :  j 

data  dt,dx,time  /0. 5, 0.0 1,50/ 


end  module  setup  info 

program  ternary_separate 

use  setup  info 
implicit  none 
interface 

!The  following  subroutine  is  designed  to  make  arrays  that 
Ifollow  the  Crank  Nicholson  Finite  Difference  method.  ’A'  will 
!have  four  columns,  representing  the  four  coefficients  in  the 
!  finite  difference  equation  mentioned  earlier  which  come  before 
!C(i-l,j+l),  C(i,j+ 1),  C(i+l,j+l),  and  C(*,j) 

!'nx’  is  the  maximum  position  of  x 

!The  'Cone'  array  represents  the  known  values  of  the  concentrations 
!  at  the  current  time  step 

subroutine  fill_a(Conc,a,nx) 
integer ::  nx 

real,  dimension(0:nx)  ::  Cone 
real,  dimension(0:nx,4)  ::  a 
!  real ::  dtt2,j 
end  subroutine  fill  a 

!The  following  subroutine  is  designed  to  solve  the  tridiagonal  matrix 
!that  was  formed  using  the  Crank  Nicholson  method 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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!It  uses  the  same  variables  as  in  subroutine  'til l  a' 

!  This  routine  returns  the  concentration  values  at  the  next  time  step 
subroutine  tridiag(Con,a,nx) 
integer ::  nx 

real,  dimension(0:nx)  ::  Con 
real,  dimension(0:nx,4)  ::  a 
end  subroutine  tridiag 

end  interface 


!This  parameter  defines  the  number  of  steps  in  the  x  direction 
!For  example,  if  the  profile  goes  from  -20  to  20,  then  NSEG=40 
integer,  parameter  ::  NSEG  =  40 
real,  parameter  ::  timestep  =  1000 
[real,  parameter ::  dt=0.5 


real,  dimension  (0:NSEG,4)  ::  a  [Coefficient  array 
real,  dimension  (0:NSEG, timestep) ::  C  [Initial  condition 
real,  dimension  (0:NSEG)  ::  xposa,conca,xposb,concb,conarray 
integer ::  count 
real ::  countreal 

[character  (len=l)  ::  tab  =  char(9) 


!***  Enter  data  and  info 

[This  statement  opens  the  text  files  containing  the 
[two  initial  concentration  profiles 
open  (unit=10,  file- trig. txt',  status-old') 
open  (unit=l  1,  file-trig2.txt',  status- old’) 

[This  statement  opens  the  ouput  file  so  the  new  concentration 

[profile  at  the  current  time  step  can  be  written 

open  (unit=100,  file-binaryoutputtrig.txt',  status- unknown') 

!dx=l/NSEG 
[print*,  'dx  is',dx 
alpha=dt/(dx*dx) 

!write(*,*)dx,dt,’  alpha  is',  alpha 

[This  loop  will  take  only  the  proper  number  of  x  positions  and 
[concentrations,  storing  them  into  'xposa'  and  'conca', 

[leaving  off  the  trailing  zeros 
count=-l; 
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5 


read  (10,*,END=15)  xposl,concl 
count=count+l; 
xposa(count)=xpos  1 ; 
conca(count)=conc  1 ; 
go  to  5 

15  if  (count.EQ.O)  then 

print*,  'No  data  in  file' 

else 
end  if 

!Same  for  second  profile 
count=-l; 

6  read  (1 1,*,END=16)  xpos2,conc2 
count=count+l; 
xposb(count)=xpos2; 
concb(count)=conc2; 
go  to  6 

16  if  (count.EQ.O)  then 

print*,  'No  data  in  file' 

else 
end  if 


[Define  the  initial  concentration  array  at  time=0  for  both  elements 


C(:,0)=concb;  !****change  a  to  b  to  a  to  change  elements 


!C(:,2,0)=concb; 
[print*,  C(:,2,0) 
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! Begin  the  loop  to  solve  over  the  entire  time,  incrementing 
!by  0.5 

t=0.0; 

countreal=0.0; 

dtt=time/.5;  ! Gives  how  large  the  time  column  is  in  the  array 
!  print*, dtt 


! print*,  C(:,2,0) 


!*****This  is  the  set  up  for  only  one  component****** 


do  j=0, 15,0.5 

!  if  (t>time)  exit 

countreal=countreal+l  hnonitor  j  increments 
!t=countreal*dt 

!  Store  one  column  of  3D  array  into  ID  array  'conarray' 

!  This  is  done  because  the  subroutine  'fill _ a'  only  needs  the 

!  concentrations  and  not  the  position  and  time  subscripts 
conarray=C(:,j) 

!  print*, z 

! print*,  conarray 
Sprint*,  C(:,  1,0) 

call  fill_a(conarray,a,NSEG)  [Fill  in  the  coefficients  for  the  current  time 


Sprint*,  C(l,i,0) 

! print*,  a(5,4) 

call  tridiag(conarray,a,NSEG)  lUpdate  C(i)  to  new  time  step 


!  This  loop  will  move  the  new  concentration  profile  back  into 
!the  3D  array  with  the  appropriate  position  and  time  subscripts 
do  k=0,NSEG 
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slot 


C(k,countreal)=conarray(k)  ! ****Won't  allow  non-integer  for  countreal 


end  do 


end  do 


! print*,  countreal 

[print*,  C(:,  1 5) !  *********  Test  print  to  show  concentration  profile  across  all  x  positions 
!  for  a  specific  element  at  a  specific  time  step  ********** 


[write  (100,99)  C(:,275);  [Write  the  current  concentration  profile  to  the  output  file 
!99  format  (IX,  F7.4); 


close  (unit=100)  [Close  the  output  file 

end  program  ternary  separate  [Close  the  program 


j ubroutine  F ill  ct 


subroutine  fill_a(Conc,a,nx) 

use  setup  info  [Uses  the  parameter  module 

integer ::  nx 
[real ::  dtt2,z 

real,  dimension(0:nx)  ::  Cone 
real,  dimension(0:nx,4)  ::  a 

[print*, Cone 

[Loop  to  enter  the  four  coefficients  of  the  finite 
[difference  equation 
do  i=l,nx-l 

a(i,l)  =  -alpha/2. 
a(i,2)  =  1 .  +  alpha 
a(i,3)  =  a(i,l) 

a(i,4)  =  Conc(i)*(l. -alpha)  +  (Conc(i-l)+Conc(i+l))*alpha/2. 

end  do 

[print*,  a(:,4); 

[These  next  two  sets  are  different  because  there  is  no  previous 


41 


Ipoint  for  i=0  and  there  is  no  future  point  for  i=NSEG 

!  This  is  for  i=0 
a(0,l)  =  0. 
a(0,2)  =  1 .  +  alpha 
a(0,3)  =  -alpha 

a(0,4)  =  Conc(0)*(l. -alpha)  +  Conc(l)*alpha 

!  This  is  for  i=NSEG 
a(nx,l)  =  -alpha 
a(nx,2)  =  1 .  +  alpha 
a(nx,3)  =  0. 

a(nx,4)  =  Conc(nx)*(l. -alpha)  +  Conc(nx-l)*alpha 

print*,  Cone 
!print*,a(nx,4); 

end  subroutine  fill  a 


|  matrix  H®H*»i*»i*»i*^*»i*^*»i*»i*»i*»i* 

subroutine  tridiag(Con,a,nx) 

real,  dimension(0:nx)  ::  Con 
real,  dimension(0:nx,4)  ::  a 
integer ::  nx 

real ::  denom  !  Represents  denominator  after  calculation 
integer ::  j 

!  Initial  values 

Con(0)=a(0,4)/a(0,2) 

a(0,3)=a(0,3)/a(0,2) 

!  This  loop  sets  all  ofthe  known  information  to  variables 
do  j=l,nx 

denom  =  a(j,2)  -  a(j,l)*a(j-l,3) 

Con(j)  =  (a(j,4)  -  a(j,l)*Con(j-l))/denom 
a(j,3)  =  a(j,3)/denom 

end  do 

!  This  loop  decides  the  next  concentration  at  the  next  time  step 
!a(nx,3)  is  0 
do  j=nx-l,0,-l 

Con(j)  =  Con(j)  -  a(j,3)*Con(j+l) 

!  print*, Con(2) 
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end  do 


end  subroutine  tridiag 
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NO.  OF 

COPIES  ORGANIZATION 


1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
ONLY)  DTICOCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FORT  BELVOIR  VA  22060-6218 

1  US  ARMY  RSRCH  DEV  & 
ENGRG  CMD 
SYSTEMS  OF  SYSTEMS 
INTEGRATION 
AMSRD  SS  T 
6000  6TH  ST  STE  100 
FORT  BELVOIR  VA  22060-5608 

1  INST  FOR  ADVNCD  TCHNLGY 
THE  UNIV  OF  TEXAS 
AT  AUSTIN 
3925  W  BRAKER  LN 
AUSTIN  TX  78759-5316 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
IMNE  ALC  IMS 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 


1  DIR  USARL 

AMSRD  ARL  Cl  OK  TP  (BLDG  4600) 
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COPIES  ORGANIZATION 


1  HQDA 

DAMO  FDT 

400  ARMY  PENTAGON 

WASHINGTON  DC  203 1 0-0460 

1  DIRECTOR 

US  ARMY  RSCH  LAB 
AMSRD  ARL  D 
J  MILLER 

2800  POWDER  MILL  RD 
ADELPHI MD  20783-1197 


2  US  ARMY  ARDEC 
AMSTA  AR  CCH  B 
C  MANDALA 
E  FENNELL 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSTA  AR  CCS 
PICATINNY  ARSENAL  NJ 
07806-5000 


1  HQDA  DIR  R&D  SAAL  TR 
W  MORRISON 
SUITE  9800 

2511  JEFFERSON  DAVIS  HWY 
ARLINGTON  VA  22201 

1  HQ  US  ARMY 
MATERIEL  CMD 
9301  CHAPEKRD 
FORT  BELVOIR  VA  22060-5527 

1  US  ARMY  BMDS  CMD 

ADVANCED  TECHLGY  CTR 
PO  BOX  1500 

HUNTSVILLE  AL  35807-3801 

1  OFC  OF  THE  PRODUCT  MGR 
SFAE  AR  HIP  IP 
R  DE  KLEINE 

PICATINNY  ARSENAL  NJ  07806-5000 

1  US  ARMY  ARDEC 

PROD  BASE  MODRNZTN  AGENCY 
AMSMC  PBM 
A  SIKLOSI 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 

PROD  BASE  MODRNZTN  AGENCY 
AMSTA  AR  WES  L  LAIBSON 
PICATINNY  ARSENAL  NJ 
07806-5000 

3  PM  PEO  ARMAMENTS 

TANK  MAIN  ARMAMENT  SYS 
AMCPM  TMA 
AMCPM  TMA  105 
AMCPM  TMA  AS 
H  YUEN 

PICATINNY  ARSENAL  NJ 
07806-5000 


1  US  ARMY  ARDEC 
AMSTA  AR  WE 
PICATINNY  ARSENAL  NJ 
07806-5000 

2  PM  MAS 

SFAE  AMO  MAS  SMC 
PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEM  D 
J  LUTZ 
BLDG  354 

PICATINNY  ARSENAL  NJ 
07806-5000 

3  US  ARMY  ARDEC 
AMSTA  AAR  AEE  W 
M  MEZGER 

D  WIEGAND 
PLU 

BLDG  3022 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AIL  F 
G  FERDINAND 
BLDG  1 

PICATINNY  ARSENAL  NJ 
07806-5000 

2  US  ARMY  ARDEC 
AMSTA  AR  WEE 
S  WESTLEY 

S  BERNSTEIN 
PICATINNY  ARSENAL  NJ 
07806-5000 
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1  US  ARMY  ARDEC 
AMSRD  AAR  AEE  W 
S  EINSTEIN 
BLDG  382 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
SFAE  AMO  CAS 
J  RUTKOWSKI 
BLDG  171  M 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEI  W 
B  BRODMAN 
BLDG  472 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEE  W 
P  OREILLY 
BLDG  3028 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
SFAE  AMO  CAS  R 
R  CIRINCIONE 
PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEE  W 
PHUI 
BLDG  382 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEE  W 
J  OREILLY 
BLDG  382 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  AMSTA  AR  FS 
TGORA 

PICATINNY  ARSENAL  NJ 
07806-5000 


1  US  ARMY  ARDEC 
AMSTA  AR  FS  DH 
PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 

AMSRD  AAR  AEE  W  F(D) 

R  KOPMANN 
BLDG  62N 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  ATD 
B  MACHAK 
BLDG  1 

PICATINNY  ARSENAL  NJ 
07806-5000 

1  US  ARMY  ARDEC 
AMSRD  AAR  AEM  C 
K CHUNG 
BLDG  407 

PICATINNY  ARSENAL  NJ  07806-5000 

1  DIR  BENET  WEAPONS  LAB 
AMSTA  AR  CCB  T 
S  SOPOK 

WATERVLIET  NY 
12189-4050 

1  DIR  BENET  WEAPONS  LAB 
AMSTA  AR  CCB  TA 
M  AUDINO 

WATERVLIET  NY  12189-4050 

1  DIR  BENET  WEAPONS  LAB 
AMST  AAR  CCB  D 

R  HASENBEIN 
WATERVLIET  NY 
12189-4050 

2  CDR  US  ARMY  RSRCH  OFC 
TECH  LIB 

D  MANN 
POBOX  12211 

RESEARCH  TRIANGLE  PARK  NC 
27709-2211 

1  PM  US  ARMY  TANK  AUTOMOTIVE 
CMD 

AMCPM  ABMS 
T  DEAN 

WARREN  MI  48092-2498 
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1  PM  US  ARMY  TANK  AUTOMOTIVE 
CMD 

FIGHTING  VEHICLES  SYSTEMS 
SFAE  ASM  BV 
WARREN  MI  48397-5000 

1  PM  ABRAMS  TANK  SYSTEM 
SFAE  ASM  AB 
WARREN  MI  48397-5000 

1  DIR  HQ  TRAC  RPD 
ATCDMA 

FT  MONROE  V A  23651-5143 
1  CDR 

RADFORD  ARMY 
AMMUNITION  PLANT 
SMCAR  QA  HI  LIB 
RADFORD  VA  24 1 4 1  -0298 

1  COMMANDANT 
USAFC&S 
ATSF  CN 
P  GROSS 

FT  SILL  OK  73503-5600 

4  CDR  NAVAL  RSRCH  LAB 
TECH  LIBRARY 

CODE  4410 
K  KAILASANATE 
J  BORIS 
E  ORAN 

WASHINGTON  DC  20375-5000 

1  OFFICE  OF  NAVAL  RSRCH 
CODE  473  J  GOLDWASSER 
800  N  QUINCY  ST 
ARLINGTON  VA  22217-9999 

5  COR  NSWC 
S  MITCHELL 
C  MICHIENZ1 
J  CONSAGA 
C  GOTZMER 
TECHLIB 

INDIAN  HEAD  MD  20640-5000 
1  CDR 

NAVAL  SURFACE  WARFARE  CTR 
CODE  G30 

GUNS  &  MUNITIONS  DIV 
DAHLGREN  VA  22448-5000 


1  CDR 

NAVAL  SURFACE  WARFARE  CTR 
CODE  G32 

GUNS  SYSTEMS  DIV 
DAHLGREN  VA  22448-5000 

1  CDR 
NSWC 
CODE  E23 
TECHLIB 

DAHLGREN  VA  22448-5000 

1  CDR 
NSWC 

R  HUBBARD  G33 
DAHLGREN  VA  22448-5000 

2  CDR 

NAVAL  AIR  WARFARE  CTR 
CODE  3895 
T  PARR 
R  DERR 

CHINA  LAKE  CA  93555-6001 
1  CDR 

NAVAL  AIR  WARFARE  CTR 
INFORMATION  SCIENCE  DIV 
CHINA  LAKE  CA  93555-6001 

1  WL  MNME 

ENERGETIC  MATERIALS  BR 
2306  PERIMETER  RD 
STE  9 

EGLIN  AFB  FL  32542-5910 

1  DIR  SANDIA  NATL  LABS 
M  BAER  DEPT  1512 
PO  BOX  5800 

ALBUQUERQUE  NM  87185 

1  DIR  SANDIA  NATL  LABS 
COMBUSTION  RSRCH  FACILITY 
R  CARUNG 

LIVERMORE  CA  94551-0469 

2  DIR  LLNL 
L355 

A  BUCKINGHAM 
M  FINGER 
PO  BOX  808 

LIVERMORE  CA  94550-0622 
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1  CIA 

J  BACKOFEN 
RM  4P07  NHB 
WASHINGTON  DC  20505 

2  MILLERSVILLE  UNIV 
PHYSICS  DEPT 

C  W  PRICE 
M  NOLAN 

MILLERSVILLE  PA  17551 

2  UNIV  OF  ILLINOIS 

DEPT  OF  MECH  INDUSTRY  ENGR 
H  KRIER 
R  BEDDINI 

144  MEB  1206  N  GREEN  ST 
URBANA  IL  61801-2978 

5  PENNSYLVANIA  STATE  UNIV 
DEPT  OF  MECHANICAL  ENGRG 
V  YANG 
K  KUO 
S  THYNELL 
G  SETTLES 
R  YETTER 

UNIV  PARK  PA  16802-7501 

1  ARROW  TECHLGY  ASSOC  INC 
1233  SHELBURNE  RD  D  8 
SOUTH  BURLINGTON  VT  05403 

1  AAI  CORPORATION 
D  CLEVELAND 

PO  BOX  126 

HUNT  VALLEY  MD  21030-0126 

2  ALLIANT  TECHSYSTEMS  INC 
ALLEGHENY  BALLISTICS  LAB 
W  B  WALKUP 

T  F  FARABAUGH 
PO  BOX  210 

ROCKET  CTR  WV  26726 

3  ALLIANT  TECHSYSTEMS  INC 
C  AAKHUS  MN07-LW54 

R  DOHRN  MN07-LW54 
D  KAMDAR  MN07-LW54 
5050  LINCOLN  DR 
EDINA  MN  55436 


4  ALLIANT  TECHSYSTEMS  INC 
RADFORD  ARMY  AMMO  PLANT 
D  A  WORRELL 
W  J  WORRELL 
S  RITCHIE 
K  BROWN 

RADFORD  VA  24141-0299 

3  ST  MARKS  POWDER 
GENERAL  DYNAMICS  ARM  SYS 
J  DRUMMOND 

J  HOWARD 
R  PULVER 
7121  COASTAL  HWY 
CRAWFORD VILLE  FL  32327 

1  GENERAL  DYNAMICS  ARM  SYS 
J  TALLEY  RM  1305 
LAKESIDE  AVE 
BURLINGTON  VT  05401 

1  PRIMEX 

BADGER  ARMY  AMMO  PLANT 
F  E  WOLF 

BARABOO  WI  53913 

4  PRIMEX 

E  J  KIRSCHKE 
A  F  GONZALEZ 
J  DRUMMOND 
D  W  WORTHINGTON 
PO  BOX  222 

SAINT  MARKS  FL  32355-0222 

2  PRIMEX 
NHYLTON  J  BUZZETT 
10101  9TH  ST  NORTH 

ST  PETERSBURG  FL  33716 

1  PAUL  GOUGH  ASSOC  INC 
P  S  GOUGH 

1048  SOUTH  ST 
PORTSMOUTH  NH  03801-5423 

2  VERITAY  TECHGY  INC 
R  SALIZONI 

J  BARNES 

4845  MILLERSPORT  HWY 
EAST  AMHERST  NY  14501-0305 
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1  PRIMEX 
E  STEINER 
DIR  LARGE  CAL  R&D 
PO  BOX  127 
RED  LION  PA  17356 

1  SRI  INTERNATIONAL 
TECH  LIB 

PROPULSION  SCIENCES  DIV 
333  RAVEN  WOOD  AVE 
MENLO  PARK  CA  94025-3493 

ABERDEEN  PROVING  GROUND 


1  CDR  USAATC 
CSTE  DTC  AT  SL 
R  HENDRICKSEN 
APG  MD  21005 

31  DIRUSARL 

AMSRD  ARL  WM 
B  RINGERS 
AMSRD  ARL  WM  BC 
M BUNDY 
J  GARNER 
P  PLOST1NS 
P  WEINACHT 
AMSRD  ARL  WM  BD 
W  R  ANDERSON 
R  A  BEYER 
A  L  BRANT 
S  W  BUNTE 
C  F  CHABALOWSKI 
T  P  COFFEE 
J  COLBURN 
P  J  CONROY 
B  E  FORCH 
B  E  HOMAN 
S  L  HOWARD 
P  J  KASTE 
A  J  KOTLAR 
C  LEVERITT 
K  L  MCNESBY 
M  MCQUAID 
M  S  MILLER 
A  W  MIZIOLEK 
J  B  MORRIS 
J  A  NEWBERRY 
M  J  NUSCA 

R  A  PESCE-RODRIGUEZ 
G  P  REEVES 
B  M  RICE 
R  C  SAUSA 
A  W  WILLIAMS 
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Intentionally  left  blank. 
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